Abstract 



In this paper, I present a new solution method for sparse regression using LO regularization. The model 
introduces a sparseness mechanism in the likelihood, instead of in the prior, as is done in the spike and slab 
model. The posterior probability is computed in the variational approximation. The variational parameters 
appear in the approximate model in a way that is similar to Breiman's Garrote model. I refer to this method 
as the variational Garrote (VG). The VG is compared numerically with the Lasso method and with ridge 
regression. Numerical results on synthetic data show that the VG yields more accurate predictions and more 
accurately reconstructs the true model than the other methods. The naive implementation of the VG scales 
cubic with the number of features. By introducing Lagrange multipliers we obtain a dual formulation of the 
problem that scales cubic in the number of samples, but close to linear in the number of features. 
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1 Introduction 

One of the most common problems in statistics is linear regression. Given p samples of 7i-dimensional input 
data x'^ji — 1, . . . ,n and 1-dimensional output data y^, with /i = 1, . . . ,p, find weights Wi, wq that best describe 
the relation 

n 

for all /i. is zero-mean noise with inverse variance /3. 

The ordinary least square (OLS) solution is given by w = X^^^ ^-nd wq — y — WiXi, where x is the input 
covariance matrix b is the vector of input-output covariances and Xi, y are the mean values. There are several 
problems with the OLS approach. When p is small, it typically has a low prediction accuracy due to overfitting. 
In particular, when p < n, x is not of maximal rank and so its inverse is not defined. In addition, the OLS 
solution is not sparse: it will find a solution Wi ^ for all i. Therefore, the interpretation of the OLS solution 
is often dificult. 

These problems are well-known, and there exist a number of approaches to overcome these problems. The 
simplest approach is called ridge regression. It adds a regularization term ^A^jiuf with A > to the OLS 
criterion. This has the effect that the input covariance matrix x g^ts replaced by x -I- A/ which is of maximal 
rank for all p. One optimizes A by cross validation. Ridge regression improves the prediction accuracy but not 
the interpretability of the solution. 

Another approach is Lasso jj. It solves the OLS problem under the linear constraint < t. This 

problem is equivalent to adding an Li regularization term A^^ \'Wi\ to the OLS criterion. The optimization 
of the quadratic error under linear constraints is called a quadratic program. There exist efficient methods to 
solve this quadratic programming problem. See [2j for a recent account. Again, A or t are found through cross 
validation. The advantage of the Li regularization is that the solution tends to be sparse. This improves both 
the prediction accuracy and the interpretability of the solution. 

The Li or L2 regularization terms are known as shrinkage priors because their effect is to shrink the size 
of Wi. The idea of shrinkage prior has been generalized by [3] to the form A^]; Iw^il' with q > and q = 1,2 
corresponding to the Lasso and ridge case, respectively. Better solutions can be obtained for q < 1, however 
the resulting optimization problem is no longer convex and therefore more difficult. 

An alternative Bayesian approach to obtain a sparse solution was proposed by [4]. They introduce n 
variational selector variables Si such that the prior distribution over Wi is a mixture of a narrow (spike) and 
wide (slab) Gaussian distribution, both centered on zero. The posterior distribution over Si indicates whether 
the solution Wi is non-zero or zero. The computation of the posterior requires the use of Gibbs sampling or any 
other MCMC method. 

Here, we propose a model where we introduce sparseness in the likelihood instead of in the prior: 

n 

y^^^Y.W^S.^x'^^+e (1) 
i=l 

with Si — 0, 1. The idea of this model is that the bits Si will identify the predictive inputs i. When we optimize 
the posterior distribution with respect to Si or using Bayesian integration, one finds the optimal subset of 

^We assume from here on without loss of generaUty that ^ X/^— i ~ p S''— i V'^ ~ ^ 
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relevant features. Since the number of subsets is exponential in n one could use MCMC sampling. Here, instead 
we propose to approximate the posterior distribution over s and w using a variational approximation. 

Model Eq.[T]has some similarity with Breiman's non- negative Garrote method [S]. The non- negative Garrote 
method assumes the same model as in Eq. [T] but with < < 1 instead of binary. It computes Wi using OLS 
and then finds Si by minimizing 

/ „ \ 2 

I — ^ x'^WiSi J subject to Si > ^ Si < t 

Because of this similarity, we refer to our method as the variational Garrote (VG). As we will see, the VG does 
not require the OLS solution and can find good solutions also when p < n. 

Using a Bayesian description, and denoting the data by D : {x^,y^},^ ~ 1, . . . ,p, the likelihood term is 
given by 

p{y\x,s,w,(3) = \/^exp f-^ _ ^ ^ 

p{D\s,w,l3) = ]^p(y''|f'',s,u;,/3) = f^^") exp ( \ ^ s^s■i'w^WjX^J - 2'^w^sA + cr^ \ U2) 




With 6, = iE,<^/^^J = |E,iyn',x^J = ^E^<<- 

We should also specify prior distributions over s, w, /?. For concreteness, we assume that the prior over s is 
factorized over the individual Si, each with identical prior probability: 

pm = f[p{s.h) pi^^h) = {^^) (3) 

with 7 given which specifies the sparsity of the solution. We denote by p{w, /?) the prior over the inverse noise 
variance /3 and the feature weights w. We will leave this prior unspecified since it can be chosen arbitrary. 
The posterior becomes 

Pis, mi) - "^"'^^"^;]]];^)'''''"'^^ (4) 

Computing the MAP estimate or computing statistics from the posterior is complex due to the discrete na- 
ture or s. We propose to compute a variational approximation to the marginal posterior piw, f3\D,^) — 
E^p(s, w, /3|I?, 7) and computing the MAP solution with respect to w, /3. Since p(D\j) does not depend on 
w, 13 we can ignore it. 

Note, that the model as specified by the likelihood Eq. [2]and the prior Eq. [3]is not equivalent to the spike 
and slab model The spike and slab model has a quadratic likelihood p{D\w,/3) that does not depend on 
selector variables. Instead, it has a prior p(iy|s) that depends on the s of the form 



/ 1 " 1 " 

\ 1=1 i=l / 

with /3± the spike and slab precisions. Thus, the spike and slab model only contains an interaction between Wi 
and Si of the form Siwf, whereas the likelihood Eq. [2] contains in addition interactions of the form WiSi. 



2 The variational approximation 

The posterior distribution Eq. 2] for given ^z;, /3 is a typical Boltzmann distribution involving terms linear and 
quadratic in s;. It is well-known that one can obtain good approximations using methods that originated in 
the statistical physics community and where Si denote binary spins. Most prominently, one can use the mean 
field or variational approximation 6 , or belief propagation . For introductions into these methods also see 
[8l[9]. Here, we will develop a solution based on the simplest possible variational approximation and leave the 
possible improvements using BP or structured mean field approximations to the future. 
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We approximate the sum by the variational bound using Jensens inequality. 



log^p(sl7M/?|s,«J,/3) >-^g(g)log , = ~Fiq,w,l3) (5) 

q{s) is called the variational approximation and can be any positive probability distribution on s and F{q, w, p) 
is called the variational free energy. The optimal q{s) is found by minimizing F[q, w, (3) with respect to q{s) so 
that the tightest bound - best approximation - is obtained. 

In order to be able to compute the variational free energy efficiently, q{s) must be a tractable probability 
distribution, such as a chain or a tree with limited tree- width |10) . Here we consider the simplest case where 
q(s) is a fully factorized distribution: q{s) — YVi=i Qiisi) with qi{si) = niiSi + (1 — nii){l — Si), so that q is fully 
specified by the expected values — qi{si = 1), which we collectively denote by m. The expectation values 
with respect to q can now be easily evaluated and the result is 

^2 ^ ^ ^ ( ^ miTOjWiWjXij + ^ mj(l - mi)w'jxii ^ 2 ^ miW.b^ + 

i.j i i—1 

- 7^m^ + nlog(l-|-exp(7)) (6) 

i=l 
n 

+ ^ (milogm^ + (1 - m^) log(l - TOi)) (7) 
1=1 

The first line is due to the likelihood term, the second line is due to the prior on s and the third line is the 
entropy of q{s). The approximate posterior marginal posterior is then 

p{w,f3\D,'y) (X p{w,f3)^p{s\j)p{D\s,w,f3) K p{w,P)exp{-F{rn,w,l3,j)) 

3 

We can compute the variational approximation rh for given w, /3, 7 by minimizing F with respect to m. In 
addition, p{w, /3\D^j) needs to be maximized with respect to w, /?. Note, that the variational approximation 
only depends on the likelihood term and the prior on 7, since these are the only terms that depend on s. Thus, 
the variational approximation does not depend on the particular choices for the prior p{w, For concreteness, 
we assume a flat prior p{w,l3) oc 1. We set the derivatives of F with respect rh,w,p equal to zero. This gives 
the following set of fixed point equations: 

f Pp 2 \ 

mi = (^\l + ^w^Xii\ (8) 



w 
1 



(xO Xtj ^ Xijmj + {1 - mj)xjjSij (9) 

n 

al -'^niiWibi (10) 



with (j{x) = (1 + exp(— a;))^^ and where in Eq. [10] we have used Eq. [9] Eqs. [8IIT0I provide the solution to the 
variational Garrote. They can be solved by fixed point iteration as outlined in Appendix [Cj Initialize rfi at 
random. Compute w by solving the linear system Eq. ^ and /3 from Eq. 1101 Compute a new solution for m 
from EqU 

Within the variational/MAP approximation the predictive model is given by 

y = ^miWiX.+S, (11) 
i 

with (^^) = 1/(3 and rh,w,(3 as estimated by the above procedure. 

Let us pause to make some observations about this solution. One might naively expect that the variational 
approximation would simply consist of replacing WiSi in Eq. [T]by its variational expectation Wiirii. If this were 
the case, rh would disappear entirely from the equations and one would expect in Eq. |9] the OLS solution with 
the normal input covariance matrix x instead of the new matrix x' (note, that in the special case that = 1 
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for all x' = X s-i^d Eq.[9]does reduce to the OLS solution). Instead, to and w are both to be optimized, giving 
in general a different solution than the OLS solution 

When TOi < 1, x' differs from x by adding a positive diagonal to it. This is similar to the mechanism of ridge 
regression, but with the important difference that the diagonal term depends on i and is dynamically adjusted 
depending on the solution for to. Thus, the sparsity prior together with variational approximation provides a 
mechanism that solves the rank problem. When all to^ < 1, x' is of maximal rank. Roughly speaking if x has 
rank p < n, x' can be still of rank n when no more than p of the to, = 1, the remaining n — p oi the to, < 1 
making up for the rank deficiency. Note, that the size of rrii (and thus the rank of x') is controlled by 7 through 
Eq.l 

2.1 Orthogonal and univariate case 

In order to obtain further insight in the solution, consider the case in which the inputs are uncorrelated: 
Xij — Sij. Then Xij — 5ij and Eq. [HI gives the solution Wi = h. Eg s. 151 and [TUl become 

1//3 

The term bfrrii is the explained variance and is subtracted from the total output variance to give an estimate 
of the noise variance 1//3. When rrii = 1, the explained variance is maximal. The Garrote solution with 
< TO, < 1 has reduced explained variance with (hopefully) a better prediction accuracy and interpretability. 
In the 1-dimensional case these equations become 

alii -mp) (13) 

with p = b'^/<Jy the squared correlation coefficient. 

In Eq. 1121 we have eliminated /3 and we must find a solution for to for this non-linear equation. We see 
that it depends on the input-output correlation p, the number of samples p and the sparsity 7. For p =^ 100, 
the solution for different p,7 is illustrated in figure [T] (see appendix Eq. 1121 has one or two solutions for to, 
depending on the values of 7, The two solutions correspond to two local minima of the free energy F , 
which in the univariate case is given by 

p 

= — log(l — pm) — 7m + const. 

where we have omitted terms independent of m. The best variational solution is given by the solution with the 
lowest free energy, indicated by the solid lines in fig.[TJ'ight. 

It is interesting to compare the univariate solution with ridge regression. Lasso or Breiman's Garrote, which 
was previously done for the latter three methods in [T]. Suppose that data are generated from the model 
y = wx -\- £, with <^^^) = ^a;^) = 1. We compare the solutions as a function of w. The OLS solution is 
approximately given by Wois ~ (xy) — w, where we ignore the statistical deviations of order 1/p due to the finite 
data set size. Similarly, the ridge regression solution is given by Wridgc ~ Aw, with < A < 1 depending on the 
ridge prior. The Lasso solution (for non-negative w) is given by luiasso = (w — 7)"*" [T], with 7 depending on 
the Li constraint. Breiman's Garrote solution is given by uigarroto = (1 ~ [T], with 7 depending on the 

Li constraint. The solutions are given in fig.[H These methods are quantitively different. The ridge regression 
solution is off by a constant multiplicative factor. The Lasso solution is zero for small w and for larger w gives 
a solution that is shifted downwards by a constant factor. Breiman's Garrote is identical to the Lasso for small 
w and shrinks less for larger w. The VG gives an almost ideal behavior and can be interpreted as a soft version 
of variable selection: For small w the solution is close to zero and the variable is ignored, and above a threshold 
it is identical to the OLS solution. 

^The technical reason that this does not occur is that in the computation of the expectation with respect to the distribution q 
that results in Eq. [7]one has {siSj) = mivrij for i ^ j, but ^s?^ = {si) = mi. 
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Figure 1: Left: Phase plot p,7 for p = 100 giving the different solutions for m. Solid and dashed lines are from 
Eq. [24] in the range of p values where two solutions for m exist. Dotted line is solution for 7 when m = 1/2, 
to indicate the transition from m w to m ~ 1. In the lower left corner, the unique solution m w is found. 
In the top right corner, the unique solution m « 1 is found. Between the dashed and the solid line, the two 
solutions TO « and to ~ 1 co-exist. Right: Example solutions for to versus p for 7 = — 10,p = 100 (top) and 
for 7 = -40, p = 100 (bottom). 
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Figure 2: Univariate solution for different regression methods. All methods yield a shrinked solution (deviation 
from diagonal line). Variational Garrote (VG) with 7 = —10 and p = 100. Ridge regression with A = 0.5. 
Garrote with 7 = 1/4. Lasso with 7 = 0.5. 
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2.2 A dual formulation 



The solution of the system of Eqs.lBlfTUlbv fixed point iteration requires the repeated solution of the n dimensional 
linear system x'w = b. When n > p, we can obtain a more efheient method using a dual formulation. 
We define new variables = rriiWix'^ and add Lagrange multipliers A'^: 

fj. i 
n n 

~ -f^rrii H- nlog(l +exp(7)) -f ^ (m^logmi + (1 - mi)log(l - m,^)) 

fi i 

We show in the Appendix [B] that the solution is found by solving the system of equations: 

p 

(15) 

(16) 

(17) 
(18) 









1 


= -E 




A^ 


= f^r 








-E 




1 


1 



— r^E^"^' (19) 



together with Eq. |S1 

2.3 Cross validation 

One remaining issue is how to choose the level of sparsity 7. If we have a prior belief about 7 we could 
optimize its value together with w, /3. Here, instead, we propose compute the VG solution for fixed 7 and 
choose its optimal value through cross validation on independent data. This has the advantage that our result 
is independent of our (possibly incorrect) prior belief. 

But another important advantage is computational. When we increase 7 from a negative value to zero in 
small steps, we obtain a sequence of solutions with decreasing sparseness. These solutions will better fit the data 
and as a result /3 increases with 7. Thus, increasing 7 implements an annealing approach where we sequentially 
obtain solutions at lower noise levels. We found empirically that this approach is effective to reduce the problem 
of local minima. To further deal with the effect of hysteresis we recompute the solution from 7max down to 

The minimal value of 7 is chosen as the largest value such that rtii = e, with e small. We find from Eqs.lSl fTOl 
that 

7mi„ = -^+a-l(6)+0(6) (20) 

We heuristically set the maximal value of 7 as well as the step size. The VG algorithm is summarized in 
Appendix 



3 Numerical examples 

In the following examples, we compare the Lasso, ridge regression and the VG. For each example, we generate 
a training set, a validation set and a test set. Inputs are generated from a mean zero multi-variate Gaussian 
distribution with specified covariance structure. We generate outputs = WiX^ + d^'^ with dif- G Af{0, a). 
For each value of the hyper parameters (7 in the case of VG, A in the case of ridge regression and Lasso), we 
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Figure 3: Top left (a): Minimal variational free energy versus 7. The two curves correspond to warm start 
solution from small to large 7 ('forward') and fr-om large to small 7 ('backward') (see also Appendix ICj) . Top 
right (b): Training and validation error versus 7. The optimal 7 minimizes the validation error, is defined as 
Bottom left (c): Solution vi = miwi and maxi=2:n I'TiiWij. The correct solution is found in the range 7 = —14.8 
to 7 « —10. Bottom right (d): Optimal solution Vi = WiUii versus i. 

optimize the model parameters on the training set. For the Lasso, we used the method described in 2 We 
optimize the hyper parameters that minimize the quadratic error on the validation set. Finally, we report the 
mean squared error on the training set, validation set and test set as well as the mean squared error and mean 
absolute error in the estimated parameters. For the VG, we define the solution as Vi = m.iWi. 

In the first example (Example 1), we take independent inputs G A/'(0, 1) and a teacher weight vector with 
only one non-zero entry: w = (1, 0, . . . , 0), n = 100 and a = 1. The training set size p = 50, validation set size 
Pa, — 50 and test set size pt = 400. We choose e = 0.001 in Eq. [20l 7max = 0.027min, A7 = — 0.027„iin- For 
details of the algorithm see AppendixO Results for a single run of the VG are shown in fig. |3l In fig. |3K, we plot 
the minimal variational free energy F versus 7 for both the forward and backward run. Note, the hysteresis 
effect due to the local minima. For each 7, we use the solution with the lowest F. In fig. we plot the 
training error and validation error versus 7. The optimal 7 — —14.8 is denoted by a star and the corresponding 
a — l/v^ = 1.0845. In fig. [3j;, we plot the non-zero component vi = miWi and the maximum absolute value 
of the remaining components versus 7. Note the robustness of the VG solution in the sense of the large range 
of 7 values for which the correct solution is found. In fig. we plot the optimal solution Vi = niiWi versus i. 

In fig. m we show the Lasso (top row) and ridge regression (bottom row) results for the same data set. The 
optimal value for A minimizes the validation error (star). In fig.|ll3,c we see that the Lasso selects a number of 
incorrect features as well. Fig. |4)d also shows that the Lasso solution with a larger A in the range 0.45 < A < 0.95 
could select the single correct feature, but would then estimate w too small due to the large shrinkage effect. 
Ridge regression gives very bad results. The non-zero feature is too small and the remaining features have large 
values. Note from fig. 0^, that ridge regression yields a non-sparse solution for all values of A. 

In tabled] we compare the performance of the VG, Lasso and ridge regression on 10 random instances. We 
see that the VG significantly outperforms the Lasso method and ridge regression both in terms of prediction 
error, the accuracy of the estimation of the parameters and the number of non-zero parameters. 

In the second example (Example 2), we consider the effect of correlations in the input distribution. Following 
[I] we generate input data from a multi-variate Gaussian distribution with covariance matrix Xij — (5'*^-'', with 
S = 0.5. In addition, we choose multiple features non-zero: Wi = l,i = 1, 2, 5, 10, 50 and all other Wi — 0. We 
use n = 100,(7 = 1 and p/pv/pt = 50/50/400. In table [5] we compare the performance of the VG, Lasso and 
ridge regression on 10 random instances. We see that the VG significantly outperforms the Lasso method and 
ridge regression both in terms of prediction error and accuracy of the estimation of the parameters. 

In fig. [S]we show the accuracy of the VG method and the Lasso method as a function of the noise a. We 
generate data according to Example 2 and vary in the range le — 4 to 10. Fig. [5] shows the error ||(5w||2 as 

^Matlab software from |http: //«w«-stat ■ Stanford. edu/-tlbs/gliimet-matlab/| 
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Figure 4: Regression solution for Lasso (top row) and ridge regression (bottom row) for same data set as in 
fig. 131 Left column: training and validation errors versus A. Middle column: Solution for the non-zero feature 
wi and the zero-features maxi=2:n \wi\. Right column: Optimal Lasso and ridge regression solution Wi versus i. 
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Val 
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# non-zero 


\\5wh 




Ridge 


0.62 ±0.43 


1.79 ±0.31 


1.77±0.15 




3.91 ±0.92 


0.79 ±0.07 


Lasso 


0.72 ±0.25 


1.09 ±0.25 


1.20 ±0.21 


10.5 ±6.5 


1.01 ±0.60 


0.19±0.11 


VG 


0.85 ±0.25 


0.97 ±0.20 


1.03 ± 0.12 


1.4 ± 0.8 


0.43 ±0.40 


0.04 ±0.05 


True 


0.93 ±0.14 


0.87 ±0.20 


0.98 ±0.04 


1 









Table 1: Results for Example 1. Train is MSE on the training set. Val is MSE on the validation set. Test is 
MSE on the test set. # non-zero is number of non-zero elements in the solution for Lasso and is X]r=i("^« ^ ^■^) 
for VG. ll'Ju'lli = J27=i l^j ~ ""^il ^i^d ||(5w||2 = X]r=i(''^« ~ ^«)^ with w the estimated solution for Ridge and 
Lasso and w — v tor VG. 





Train 


Val 


Test 


# non-zero 


\\6w\\i 


\\Sw\\2 


Ridge 


0.29 ±0.31 


3.52 ±0.69 


3.48 ±0.35 




11.3 ±0.69 


2.79 ±0.13 


Lasso 


0.78 ±0.39 


1.40 ±0.33 


1.47 ±0.21 


16.4 ± 6.0 


2.00 ±0.64 


0.51 ±0.19 


VG 


0.83 ±0.26 


1.17±0.20 


1.21 ±0.17 


4.9 ±0.3 


0.96 ±0.48 


0.30 ±0.36 


True 


0.93 ±0.14 


0.87 ±0.20 


0.98 ±0.04 
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Table 2: Results for Example 2. Train is MSE on the training set. Val is MSE on the validation set. Test is 
MSE on the test set. ^ non-zero is number of non-zero elements in the solution for Lasso and is Yl^=ii''^i ^ ^■^) 
for VG. ||(5it;||i = X]"=i l""^* ~ ""^il ^^-^ ||i5w||2 — ^ with w the estimated solution for Ridge and 

Lasso and w — v for VG. 
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Figure 5: Accuracy of VG and Lasso in terms of the error ||(5u'||2 as a function of the noise. Data according to 
Example 2 with 6 = 0.5.dd AU results are averages over 10 runs. The variance in these results is of order 50 







Example la 






Example lb 








\\dw\\2 


max(|w3|) 


\\dw\\i 


\\dw\\2 


max(|w3|) 


Ridge 


0.6366 ±0.1811 


0.1504 ± 0.0787 


0.48 


0.1905 ±0.1113 


0.0185 ±0.0206 


0.27 


Lasso 


0.1916 ±0.1363 


0.0199 ± 0.0260 


0.30 


0.0628 ±0.0364 


0.0030 ±0.0032 





Garrote 


0.0491 ± 0.0257 


0.0020 ±0.0018 


2e- 14 


0.0491 ±0.0257 


0.0020 ±0.0018 


2e- 14 



Table 3: Accuracy of Ridge, Lasso and Garrote for Example la,b from pT. P ^ Pv ~ 1000. Parameters A (Ridge 
and Lasso) and 7 (Garrote) optimized through cross validation. ||rfu'||i^2 a-s before, maxdwsl) is maximum over 
100 trials of the absolute value of w^. Example la is inconsistent for Lasso and yields much larger errors than 
the Garrote. Example lb is consistent and the quality of the Lasso and Garrote are similar. Ridge regression 
is bad for both examples. 

defined above. We distinguish three noise domains: for very large noise both VG and LASSO produce errors 
of 0(1) and fail to find the predictive features. For intermediate and low noise levels, VG is approximately 1 
order of magnitude more accurate than Lasso. In the limit of zero noise, the error of VG keeps on decreasing 
whereas the Lasso error saturates to a constant value. 

It is well-known that the Lasso method may yield inconsistent results when input variables are correlation. 
[llj derive a necessary and sufficient condition for consistency. In addition, they give a number of examples 
where Lasso gives inconsistent results. Their simplest example has three input variables, xi,X2,X3. xi,X2,£,,e 
are independent and Normal distributed random variables, X3 = 2/3xi + 2/3x2 + £, and y = ^i^i + e, 

p = 1000. When w = (2,3,0) this example violates the consistency condition. The Lasso and VG solutions 
for different values of A and 7, respectively are shown in fig. |6l The average results over 100 instances of these 
problems is shown in table [3] 

In fig. [7] we show how the Lasso, the ridge regression and the VG behave as a function of the number of 
samples. The VG is significantly more accurate than the Lasso, for almost all values of p. This is evident from 
the ||(5if II2 and validation error. This is also evident by comparing the Lq and Li norms of the solutions found. 
The Lasso solution find a solution with the same (correct) Li norm as the VG, but spreads this solution over 
many more features, as is witnessed by the Lq norm. 

In fig. [8] we show how the Lasso, the ridge regression and the VG behave as a function of the number of 
samples. For the VG we use the dual method described in section [221 The VG has constant quality in terms of 
the error ||5w||2 and close to optimal norms Lq = Li = 5. The computation time of the VG scales approximately 
linear with n. The Lasso is significantly faster, but its quality deteriorates with n. 
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Figure 6: Top left) Lasso solution versus A is called inconsistent because it does not contain a A for which the 
correct sparsity w(S) = is obtained. Top right) Garrote solution for v versus 7 contains large range of 7 
for which the correct solution is obtained. Bottom left) Garrote solution for m. Curves for 7711^2 are identical. 
Bottom right) Garrote solution for w. 
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Figure 7: Performance of different methods as a function of the number of samples p. Data are gener- 
ated as in Example 2. n = 500 features and (3 = 1,(5 — 0.5. The number of samples is varied from 
p — [10,20,50,100,250,500]. Optimal hyper parameters found through cross validation. Cross validation 
set of size p„ = 50. Results are average over 10 runs. Top left; Quadratic error ||(5w||2 — X]"=i(^i ~ WiY' . Top 
right: Lq norm of the solution found. For the VG the number of non-zero components is defined as those i 
for which > 0.5. Bottom left: Li norm of the solution found. For the VG the Li norm of the vector with 
components Wivrii is used. Bottom right: validation error. 
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Figure 8: Performance of different methods as a function of the number of features n. Data are generated as 
in Example 2. p — 100, /3 = 2,5 = 0. Optimal hyper parameters are found through cross vahdation. Cross 
vahdation set of size Pu = 100. The number of features is varied in steps of 100 from n — 100 to n = 1000, 
and n = 2000 and n = 5000. For each value averages over 10 runs are plotted. Top left: Quadratic error 
ll'^'^^IU = X]r=i(^« ~ Wi)'^. Top right: Lq norm of the solution found. For the VG the number of non-zero 
components is defined as those i for which nii > 0.5. Bottom left: Li norm of the solution found. For the VG 
the Li norm of the vector with components Wirrii is used. Bottom right: Matlab CPU time in seconds for VG 
and C implementation CPU time in seconds for Lasso. 
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4 Discussion 



In this paper, I have introduced a new model for sparse regression by introducing binary switch variables s that 
muhiply the feature weights. When using a Bayesian framework one obtains a posterior distribution that is 
intractable due to the binary nature of s. In this paper we have proposed to approximate the variables s using a 
variational approximation. We have demonstrated that this approach is efficient and accurate and yields results 
that significantly outperform the Lasso approach. 

There are many points that need further investigation. The variational approximation can by replaced by 
any other approximate inference method. It would be interesting to see whether the results could be further 
improved by using Belief Propagation, the Cluster Variation Method or Monte Carlo Sampling. 

I have not explored the use of different priors on w or on (3. In addition, a prior could be imposed on 7. It 
is likely that for particular problems the use of a suitable prior could further improve the results. 

One of the important advantages of the Lasso over the VG is that the Lasso problem is convex whereas the 
VG may have local minima, as is illustrated in fig. [TJ This is clearly a disadvantage of the VG, but is shared 
by all other methods that use Lq regularization. We have not much explored the issue of local minima. For 
the examples that we presented the annealing approach that results from increasing 7, followed by a 'heating' 
phase to detect hysteresis, seems to work well. We have noted that for problems with a large number of features 
{0{n)), small p and high noise this approach may fail. The local minima problem depends on the approximation, 
and may be more or less severe when using other approximation schemes. 

For efficient implementation one may use the fact that the data matrix is sparse. For instance in genetic 
regression problems this may be the case, when Xi denotes single nucleotide polymorphisms (SNPs) that have 
values Xi = 0,1,2. For each i one may code them such that the value = is the most frequent one. In 
that case, substracting the mean from the data such that = will yield x^ non-sparse. The model as 

presented here should then be slightly generalized to estimate a constant bias as well. 

In principle, the model Eq. [l]when used in a full Bayesian setting is non-linear: The predictive distribution 



is a mixture of 2" Gaussians. Only in the variational/MAP approximation does this model reduce to the linear 
model Eq. [TlJ It would be of interest to explore the Bayesian extension of the VG. However note, that in our 
experiments we find that rm « or « 1 most of the time and encounter very few intermediate or multi- modal 
solutions. This suggest that the Bayesian integration, at least as far as the binary variables is concerned may 
be of limited additional value for the examples considered in the paper. 

We have seen, that the performance of the VG is excellent also in the zero noise limit. In this limit, the 
regression problem reduces to a compressed sensing problem [121 113j . The performance of compressed sensing 
with Lq sparseness penalty was analyzed theoretically in ^T^, showing the superiority of the Li penalty in 
comparison to the L2 penalty and suggesting the optimality of the Lq penalty. Our numerical result are in 
agreement with this finding. 
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A Phase plot computation for the univariate case 

/(m) in Eq. [T2]is an increasing function of m and crosses the line m either 1 or three times, depending on the values of 
p and 7 (see fig. [9)l . 

Thus, there are regions of parameter space "/,p,p where the univariate solution is unique and others for which there 
are two solutions. The transition between these two regions is when f'{m) — 1 and /(m) = m. These two equations 
imply 

(1 + ipp')m2 + {-2p - ^pp^)m +1 = (22) 

This quadratic equation in m has either zero, one or two solutions, corresponding to no touching, touching once and 
touching twice, respectively. Denote a — 1 + ^pp^ ,h — 2p + \pp^ . The critical value for p,p is when Eq. [22] has one 
solution for m, which occurs when 

D = 6^ - 4a = ipV - 2p^(l - p)p + 4(p^ - 1) = 

This is a quadratic equation in p and we solve for p in terms of p: 

P-^VW(VW±V2) 
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Since p < 1 there is a positive and a negative solution for p. The latter we disregard. 

:^/l~p(yi~p+^) (23) 



4 

p2 



p* is a decreasing function of p. 

For p > p* , D > and we find two solutions for m, which we denote by mi, 2 = Note, that the solutions in 

these critical points only depend on p,p. For each of these solutions we must find 7 such that /(m) — m, which is given 
by 

7, = log^-f^^ ..1,2 (24) 
1 — mi 2 1 — pnii 

For p < p*, 7? < and we find no solutions for m. In this case the conditions f'{m) — 1 and /(m) = m cannot be jointly 
satisfied. 

We illustrate the phase plot p, 7 for p — 100 in fig. [1] 

B Derivation of dual fixed point equations 

We compute the derivatives of Eq. 1141 



dF 

dwi 

dF 



mi /3p(l - mi)xiiWi 



= /3(z'^ - y-) + A" 



= -2^+22^^^ +-}^m.(l-mO».Xi. 
= ^(1 - 2mi)wiXii - 7 + o-"'^(mi) - 



OF 



By setting = = we obtain Eq. llQI and — — \ . Setting the remaining derivatives to zero, and eliminating 



Wi and 2:'' we obtain Eq. [S] and 



/3 = iVA^A^A^, (25) 

p ^ ^ 

/Sy" = E^M^A" (26) 

V 

with given by Eq. (TH] For given yl/j^, let ij denote the solution of Ay — y. Then it is easy to verify that Eqs. I16I17I 
solve the system of Eqs. I25I26I 
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C The Variational Garrote algorithm 

Here we give a summary of the VG algorithm. 

1. Preprocess data x'^ ,fj,= 1, . . . , p such that — — 0. 

2. Compute b. = i^^<y" 

3. For n <p, compute Xij = ^ Z^p ^i^j- 

4. Choose e and compute jmin from Eq. [20l Choose A7 and 7max. 

5. Initialize all = 0. 

6. For 7 = 7min : A7 : 7max do 

(a) Set rj = l. 

(b) While not converged: 

i. Compute w, (3 from Eqs. [QllTOl (n < p) or Eqs. I15I19I (n > p) . 

ii. Compute m using a smoothed version of Eq. |8] = (1 — 7;)mi + r/(7(. . .). 

iii. When maxi — m^j > 0.1 divide by two. 

iv. m = m'. 

(c) Store solution as wi, mi, /3i I7. 

(d) Compute -^1(7) = F(it)i, mi, /3i I7) from Eq. [3 

7. For 7 = 7niax : -A7 : 7niin do 

(a) As 6 (a) 

(b) As 6 (b) 

(c) Store solution as i(J2, m2, /32I7. 

(d) Compute -^2(7) = -F(w)2, m2, /32I7) from Eq. [3 

8. For 7 = 7niin : A7 : 7max 

(a) choose solution w, m, /3|7 that has minimal -Fi,2(7). 

(b) compute cross validation error, defined as the quadratic error on validation data using Eq. [11] 

9. Output solution w, rn, /?, 7 with minimal cross validation error. 
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